home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_09_05 / 9n05050a < prev    next >
Encoding:
Text File  |  1991-02-05  |  2.6 KB  |  157 lines
  1. /*  fast Fourier transform (FFT)
  2.  
  3.  
  4. from Handbook of C tools for Scientists and Engineers by L. Baker
  5.  
  6. DEPENDENCIES:
  7.  
  8. header file complex.h required
  9.  
  10. */
  11. #include <stdio.h>
  12. #include "complex.h"
  13. #define twopi 6.283185307
  14.  
  15. int iterp;/* global to return count used*/
  16.  
  17. #define max(a,b) (((a)>(b))? (a): (b))
  18. #define abs(x) ((x)?  (x):-(x))
  19. #define DOFOR(i,to) for(i=0;i<to;i++)
  20.  
  21. main(argc,argv) int argc;char **argv;
  22. {int i,ii,nh,n=16;
  23. double invn,exp(),dt=.25,omega,realpt,imagpt;
  24. struct complex w[32],wi[32], data[32];
  25. nh=n>>1;
  26. DOFOR(i,n)
  27.     {
  28.     CMPLX(data[i], exp(-dt*(i)),0.);
  29.     };
  30. /* caveat*/ data[0].x=.5;/*not 1. see text*/
  31. printf(" before transform:\n");
  32. DOFOR(i,n){printc(&(data[i])) ;printf("\n");}
  33. invn=1./n;
  34. fftinit(w,wi,n);
  35. printf(" w factors:\n");
  36. DOFOR(i,n){printc(&(w[i])) ;printc(&(wi[i])) ;printf("\n");}
  37. printf(" transformed:\n");
  38. fft(data,w,n);
  39. DOFOR(i,n){printc(&(data[i])) ;printf("\n");}
  40. printf(" scaled by dt and compared to analytic answer:\n");
  41. DOFOR(i,n)
  42.     {
  43.     ii= (i-nh) ? i : n-i ;
  44.     omega= twopi*ii/(n*dt);
  45.     realpt= 1./(1.+omega*omega);
  46.     imagpt= -realpt*omega;
  47.     printf(" %f %f  %f %f\n", dt*data[i].x,dt*data[i].y,realpt,imagpt);
  48.     }
  49. /* inverse fft*/
  50. fft(data,wi,n);
  51. printf(" transformed back and scaled by 1/N:\n");
  52. DOFOR(i,n)
  53.     {
  54.     CTREAL(data[i],data[i],(invn));
  55.     printc(&(data[i]));printf("\n");
  56.     }
  57. exit(0);
  58. }
  59.  
  60.  
  61.  
  62. int bitr(k,logn) int k,logn;
  63. {
  64. int ans,j,i;
  65. ans=0;
  66. j=k;
  67. DOFOR(i,logn)
  68.     {
  69.     ans=(ans<<1)+(j&1);
  70.     j=j>>1;
  71.     }
  72. return(ans);
  73. }
  74.  
  75.  
  76. int log2(n) int n;
  77. {
  78. int i;
  79. i=-1;/* will return -1 if n<=0 */
  80. while(1)
  81.     {
  82.     if(n==0)break;
  83.     n=n>>1;
  84.     i++;
  85.     }
  86. return(i);
  87. }
  88.  
  89. printc(x) struct complex *x;
  90. {
  91. printf("%f %f",x->x,x->y);
  92. return;
  93. }
  94.  
  95. fftinit(w,wi,n) int n; struct complex w[],wi[];
  96. {
  97. int i;
  98. double realpt,imagpt,cos(),sin();
  99. double factr,angle;
  100. factr=twopi/n;
  101. DOFOR(i,n)
  102.     {angle=i*factr;
  103.      realpt=cos(angle);imagpt=sin(angle);
  104.     CMPLX(w[i],realpt,imagpt);
  105.     CMPLX(wi[i],realpt,(-imagpt));
  106.     }
  107. return;
  108. }
  109.  
  110.  
  111. fft(x,w,n) int n; struct complex w[],x[];
  112. {
  113. int n1,logn,i,j,k,l,logl,p;
  114. struct complex s,t;
  115. logn=log2(n);
  116. n1=n>>1;
  117. j=logn-1;
  118. /* transform*/
  119. k=0;
  120. DOFOR(logl,logn)
  121.     {
  122.     do{
  123.     DOFOR(i,n1)
  124.         {
  125.         p=bitr((k>>j),logn);
  126.         l=k+n1;
  127.             CONJG(s,w[p]);
  128.             CMULT(t,s,x[l]);
  129.             CSUB(x[l],x[k],t);
  130.             CADD(x[k],t,x[k]);
  131.             k++;
  132.         };/* dofor i*/
  133.     k+=n1;
  134.     }while (k<n);
  135. k=0;
  136. j--;
  137. n1=n1>>1;
  138. }
  139.  
  140.  
  141. /*reorder*/
  142. for(i=1;i<n;i++)
  143.     {
  144.     k=bitr(i,logn);
  145.     if (i>k)
  146.         {
  147.         /*exchange i,k elements*/
  148.         CLET(s,x[i]);
  149.         CLET(x[i],x[k]);
  150.         CLET(x[k],s);        
  151.         }
  152.     
  153.     };
  154.  
  155. return;
  156. }
  157.